Transcriptomics of ivermectin response in Caenorhabditis elegans: Integrating abamectin quantitative trait loci and comparison to the Ivermectin-exposed DA1316 strain

Parasitic nematodes pose a significant threat to human and animal health, as well as cause economic losses in the agricultural sector. The use of anthelmintic drugs, such as Ivermectin (IVM), to control these parasites has led to widespread drug resistance. Identifying genetic markers of resistance in parasitic nematodes can be challenging, but the free-living nematode Caenorhabditis elegans provides a suitable model. In this study, we aimed to analyze the transcriptomes of adult C. elegans worms of the N2 strain exposed to the anthelmintic drug Ivermectin (IVM), and compare them to those of the resistant strain DA1316 and the recently identified Abamectin Quantitative Trait Loci (QTL) on chromosome V. We exposed pools of 300 adult N2 worms to IVM (10−7 and 10−8 M) for 4 hours at 20°C, extracted total RNA and sequenced it on the Illumina NovaSeq6000 platform. Differentially expressed genes (DEGs) were determined using an in-house pipeline. The DEGs were compared to genes from a previous microarray study on IVM-resistant C. elegans and Abamectin-QTL. Our results revealed 615 DEGs (183 up-regulated and 432 down-regulated genes) from diverse gene families in the N2 C. elegans strain. Of these DEGs, 31 overlapped with genes from IVM-exposed adult worms of the DA1316 strain. We identified 19 genes, including the folate transporter (folt-2) and the transmembrane transporter (T22F3.11), which exhibited an opposite expression in N2 and the DA1316 strain and were deemed potential candidates. Additionally, we compiled a list of potential candidates for further research including T-type calcium channel (cca-1), potassium chloride cotransporter (kcc-2), as well as other genes such as glutamate-gated channel (glc-1) that mapped to the Abamectin-QTL.


Introduction
Parasitic nematodes cause chronic and debilitating illnesses in humans and animals, as well as economic losses in livestock and agriculture [1][2][3]. For example, soil-transmitted helminths, infect more than one billion people and cost roughly five million disability-adjusted life years reported in 2010 [4] and plant-parasitic nematodes affecting crops cause annual economic losses of >$80 billion globally [5]. Control strategies for parasitic nematodes in domesticated or livestock animals rely on anthelmintic drugs, which cost the European ruminant livestock industry an estimated €320 million per year [6]. The three major classes of anthelmintic drugs are benzimidazoles, tetrahydropyrimidines, and macrocyclic lactones (MLs). The MLs, such as Ivermectin (IVM) and Abamectin, are the most commonly used drug class in veterinary/ human medicine and agriculture due to their high efficacy, low toxicity, and broad-spectrum nature [7]. Glutamine-gated chloride ion channels (GluCls) are the primary target of MLs in nematodes. This consequently leads to the depolarization of pharyngeal muscles and hyperpolarization of body wall muscles [8]. As a result, feeding becomes impaired, and the worm experiences flaccid muscle paralysis [9,10]. In addition, γ-aminobutyric acid (GABA) receptors have been reported as secondary targets of IVM at concentrations [11] considered therapeutically insignificant.
However, due to the widespread use of anthelmintic drugs, extensive anthelmintic drug resistance has developed in many parasitic nematode species of livestock [12]. Anthelmintic resistance has been reported in Haemonchus contortus, Cooperia spp., and Ostertagia ostertagi infecting sheep in Europe [13]. In addition, anthelmintic-resistant parasitic nematodes in sheep, cattle and goats have been reported in a number of countries, including Brazil, Argentina, the United States, Australia, New Zealand, and South Africa (reviewed [12]). Furthermore, there are growing concerns that mass drug administration (MDA) programs are selecting for anthelmintic resistance among helminths that infect humans [14,15]. In recent years, there have been emerging reports of IVM resistance in human parasite Onchocerca volvulus in Ghana [16] and Cameroon [17]. Evidently, understanding and identifying the genetic markers underlying resistance in parasitic nematodes is necessary.
However, identification of candidate genes or gene variants that drive anthelmintic resistance in parasitic nematodes can be challenging. This is partly due to the fact that parasitic nematodes have a host-dependent life cycle, are difficult to cultivate, maintain, manipulate, and study at scale, and often lack well-annotated genomes in full chromosomes. Until now, only a few parasitic nematodes, have relatively complete chromosomal-level genomes [18][19][20][21], allowing genetic and comparative genome analyses. In comparison, the free-living nematode Caenorhabditis elegans is biologically simple, has a rapid 3.5-day life cycle, genetically amenable, and molecular tools for manipulation are readily available [22]. In addition, the evolutionary relationship between C. elegans with parasitic nematodes [23], has led to the adoption of this free-living nematode as a model for parasitic nematodes [22,24] because it allows for comparative studies and translatable results [25].
Studies in C. elegans have begun to uncover IVM resistance mechanisms. For example, simultaneous mutations in glc-1, avr-14, and avr-15 result in a~4000-fold resistance to IVM in C. elegans [26]. These genes, together with glc-2, glc-3, and glc-4, code for the GluCls subunits, which are capable of forming diverse homomeric and heteromeric chloride channels [9,10]. However, the discovery of IVM-resistant populations of the small-ruminant parasite H. contortus lacking mutations in the glc-1, avr-14, and avr-15 genes [19,27] suggests the existence of alternative mechanisms of resistance to MLs. Several resistance mechanisms to MLs have been proposed, including changes in gene expression in GluCl and transport protein genes [28]. For example, as previously postulated, decreased expression of drug target genes may result in fewer drug receptors and thus reduced drug efficacy [29,30]. Gene expression changes of transporter genes of the ATP-binding cassette family, particularly Pglycoproteins (Pgps), also known as efflux pumps, have been linked to anthelmintic drug resistance [31,32]. This is because efflux pumps eliminate drugs from the cell, preventing them from reaching their target sites. Furthermore, increased expression of genes encoding drug metabolic enzymes, postulated to increase conversion of drugs into inactive metabolites, has been associated with resistance [33]. Collectively, this highlights the multigenic nature and multiple mechanisms of anthelmintic resistance. Although a previous microarray-based study [34] examined the effect of IVM on gene expression in C. elegans, the strain used was DA1316, which is highly resistant to IVM due to simultaneous mutations in three genes (glc-1, avr-14, and avr-15). As a result, the current study is premised on the lack of a previous transcriptomic study assessing the effect of IVM on gene expression in wild-type (N2) C. elegans.
In this study, we investigated gene expression in an adult C. elegans N2 strain exposed to IVM using an RNAseq transcriptomic approach in order to identify the underlying processes and genes involved in the response to IVM in a wild type strain. In addition, we compared our transcriptomic data with microarray data from the IVM-exposed C. elegans DA1316 strain [34] in order to identify genes that may be uniquely associated with IVM response in IVMresistant strains. Additionally, we investigated the recently discovered Abamectin Quantitative Trait Loci (QTL) on chromosome V in C. elegans [35] to determine any potential associations between the identified genes and the QTL.
A detailed description of the experimental setup is described in [41]. Briefly, synchronized worms at maximum reproduction stage (76 h post L1) were incubated at 20˚C in S-complete media supplemented with IVM in concentrations of 10 −7 and 10 −8 M (+ 0.025% DMSO final concentration) and 0.025% DMSO (control) for 4 h. The experiment was performed in quadruplicates (~300 worms /replicate). After exposure worms were pelleted by centrifugation, frozen in liquid nitrogen, ground with a pestle and subsequently suspended in Trizol (Invitrogen, Carlsbad, USA). Chloroform was added and the aqueous phase was advanced to NucleoSpin 1 RNA Plus Kit (Macherey Nagel, Düren, Germany) for RNA extraction. RNA quality and quantity checks were performed using the RNA ScreenTape kit on TapeStation 4150 (Agilent, Santa Clara, USA).

Library preparation and RNA sequencing
The SNP -& SEQ Platform, SciLifeLab Uppsala, Sweden, did the library preparation and sequencing. Using the TruSeq stranded mRNA library preparation kit with polyA selection (Illumina Inc, San Diego, USA); one microgram of RNA from each sample (biological replicate) was used to prepare sequencing libraries. From the libraries, clusters were made, and 150 cycles of paired-end sequencing were done in a single-end SP flowcell using NovaSeq 6000 equipment and v1.5 sequencing chemicals (Illumina Inc., San Diego, USA).
Expression quantification against the C. elegans transcriptome was performed with Salmon (v1.9.0) [47] using the entire genome (PRJNA13758, release WS283) as the decoy sequence. The transcriptome used was a concatenation of mRNA, ncRNA, and pseudogene files of PRJNA13758, release WS283 available at WormBase [48]. We have made the RNA sequencing raw files available at the ENA database, accessible through the link https://www.ebi.ac.uk/ena/ browser/view/ (accession number: PRJEB59331). Additionally, we have uploaded the pipeline for analysis and the necessary files for read processing, mapping, and quantification to the following GitHub repository: https://github.com/ruqse/N2IVM. R version 4.1.2 was utilized for all R-based analyses. Salmon output was matrix-summarized with tximport R package (v1.22.0) [49] utilizing the TxDb.Celegans.UCSC.ce11.ensGene annotation R package (v3.12.0) [50], for downstream transcript/ gene-level analysis. As a quality control measure, exploratory analysis and visualization was performed on count data.

Principal component analysis
Principal component analysis (PCA) was utilized to investigate sample variation. The variance stabilizing transformation function from the DESeq2 R package (v1.34.0) [51] was applied to count data of genes across all samples (n = 12). Since IVM exposure was performed in batches based on concentration, any batch effects were removed by removeBatchEffect function from the limma R package (v3.50.3) [52]. Using R base prcomp function, principal component analysis was performed on the variance-stabilized transformed (VST) read counts. The top two PCs that account for the most variation in the data were visualized in PCA plot for genes with non-zero read counts using the ggplot2 R package(v3.3.5) [53]. Sample clustering was determined by the enclosing ellipses implemented by the geom_mark_ellipse function in the ggforce R package (https://github.com/thomasp85/ggforce/).

Differential gene expression
Differential gene expression was determined on raw count data by DEseq function in the DESeq2 R package using the following model design,~batch + IVM concentration. The batch variable was to compensate for any biases caused by the batch effect. Differentially expressed genes (DEGs) were defined as those with an adjusted P-value < 0.05 and Log2foldchange (Log2FC) � 0.5 or � -0.5 (1.4-fold change). These cut-offs are based on recommendations by Schurch et al. [54]. Differentially expressed genes with average expression level (baseMean) of <20 across all samples were discarded based on the independent filtering results from DESeq2.
The metadata of DEGs such as concise description, gene ontology association, interacting gene, tissue expression and RNAi/allele phenotype observed were retrieved from WormBase using the SimpleMine tool taking DEGs gene IDs as input. To assess genes shared by each drug concentration, comparative Venn diagrams were constructed using gene IDs and the venn.diagram function of the VennDiagram R package [55]. The get.venn.partition function was utilized to acquire the gene IDs associated with each Venn diagram partition. The hypergeometric test was used to test the significance of the overlapping genes between concentrations using the R base function, phyper.

Over-representation analysis (ORA)
Over-representation analysis was performed on DEGs using all genes (19565) with a non-zero total read count as a background based on Gene Ontology (GO: Release 2021-12-15), Kyoto Encyclopedia of Genes and Genomes (KEGG: Release 2021-12-27) and TRANSFAC (Release 2021.3 classes: v2) databases. The analysis was carried out using the gost function from the gprofiler2 R package (v0.2.1) [56], with default parameters and C. elegans as the set organism. The output files were processed, and data visualized using ggplot2. Differentially expressed genes from enriched terms putatively involved in drug metabolism/resistance, apoptosis and transcription of regulation were retrieved, categorized, and characterized.

Comparative analysis gene expression between C. elegans strains N2 and DA1316 after IVM exposure
A comparative analysis of the obtained DEGs with those from a previous study [34] was performed. In that study, C. elegans strain (DA1316), a triple mutant of the glutamate-gated chloride channel (GluCI) subunits avr-14, avr-15, and glc-1 was exposed to IVM 10 −6 and 10 −7 M for 4 h. Differential gene expression was determined by microarray assays. To reduce ambiguity prior to comparison, gene names and IDs of DEGs from the previous study were evaluated to determine if they had been renamed, merged, or divided, using the WormBase Gene Name Sanitizer tool. Only DEGs with a fold change of at least 1.4 were included in the comparison and visualized in upsetR plot using the UpSetR R package (v1.4.0) [57]. The hypergeometric test was used to test the significance of the overlapping genes between DA1316 and N2 strain using the R base function, phyper.

Analysis of IVM-induced gene expression within Abamectin-QTL on chromosome V
A recent study described quantitative trait loci (QTL) regions on chromosome V of C. elegans after exposure to the macrocyclic drug, Abamectin [35]. Based on genome wide association mapping, two regions were detected, left arm (VL) and right arm (VR) located at nucleotides 1,757,246-4,333,001 and 15,983,112-16,599,066, respectively. Linkage mapping approach identified three regions, VL, VR and central (VC) located at 2,629,324-3,076,312, 15,933,659-16,336,743 and 6,118,360-7,374,129, respectively. In our comparative analysis, VL and VR from the genome wide association mapping and VC from the linkage mapping were used. The BSgenome.Celegans.UCSC.ce11genome object from the BSgenome.Celegans.UCSC.ce11 R package (v1.4.2) [58] was used in this analysis. All genes in the VL, VC and VR regions were extracted using the GRanges function and subset based on our DEGs using subsetByOverlaps function. Both functions are derived from GenomicRanges R package (v1.46.1) [59]. To eliminate any chromosome length biases, the distribution/ enrichment of DEGs across all chromosomes was estimated using the Fisher's Exact Test for Count Data in R. Differentially expressed genes mapping to the QTL were subsequently identified and visualized on a karyoplot from karyoploteR R package (v1.20.3) [60].

Transcriptomic variation amongst adult C. elegans worms exposed to Ivermectin
Total RNA (RIN > 8) was sequenced from C. elegans adults at maximum reproduction exposed to IVM 10 −7 and 10 −8 M yielding sequencing data of 21-60 million read-pairs/sample. The used concentrations of IVM treatment were selected according to [41]. However, due to the very low number of differentially expressed genes in IVM 10 −9 M, only data from IVM 10 −7 and 10 −8 M was used in the current study. Following quality control and evaluation with SortMeRNA and Trimmomatic, 16-42 million read-pairs/sample were obtained with a mean phred score above 33, indicating high quality reads. The read-mapping rate to the custombuilt PRJNA13758 C. elegans transcriptome (see Methods) and the C. elegans genome (PRJNA13758, release WS283) as a decoy, was 97% using Salmon.
According to principal component analysis of the VST read counts, worm pools (n = 300 worms per pool), both PC1 and PC2 explained the separation by treatment condition (Fig 1).

PLOS ONE
Transcriptomics of Ivermectin response in C. elegans: Integrating abamectin QTLs and DA1316 strain diverse set of genes or gene families including transmembrane transporters, and heat shock proteins (S3 Table).
To gain insights into the biological pathways underlying the differentially expressed genes, we conducted over-representation analysis on both the upregulated (Fig 3A & S4 Table) and downregulated genes (Fig 3B & S4 Table). We found that the downregulated DEGs exhibited more enrichment terms (42) compared to the upregulated genes (17). The ORA on the upregulated DEGs revealed enrichment of terms related to response to external stimulus, ABCtransport activity, and Factor: elt-3; motif: TCTTATCA (TF: M07154) based on GO, KEGG pathway, and TRANSFAC databases, respectively (Fig 3A & S4 Table). In contrast, ORA on the downregulated genes showed enrichment of terms related to metabolism of glycans, amino acids, and organic compounds, Factor: elt-3; motif: TCTTATCA (TF: M07154) among others, based on the same databases ( Fig 3B & S4 Table). Differentially expressed genes from ORA terms with putative involvement in drug metabolism were explored further.

Differentially expressed genes encoding drug receptors.
After exposure of N2 worms to IVM 10 −7 M, the putative drug target glc-1, an alpha subunit of the GluCI receptor, glc-1 was downregulated 1.4-fold. In addition, the putative drug target, lgc-26, a possible

PLOS ONE
nicotinic acetylcholine receptor was downregulated 1.6-fold after exposure to IVM 10 −7 M (Table 1). However, in IVM 10 −8 M concentration, differential expression of the putative drug receptors was not observed.

PLOS ONE
Transcriptomics of Ivermectin response in C. elegans: Integrating abamectin QTLs and DA1316 strain

Ivermectin-induced DEGs are mapped within the Abamectin QTL on chromosome V
Evans et al. [35] identified three major QTL (VL, VR and VC) correlated to natural variation in Abamectin response, on chromosome V in C. elegans. We decided to use this information to see if any of the DEGs induced by IVM, map in the QTL regions. There was significant (adjusted p = 0.01) enrichment for DEGs in chromosomes X and V (S6 Table). A comparison of the DEGs of IVM-exposed C. elegans N2 adult worms to the QTL revealed that 45 DEGs Table 3. Overlapping differentially expressed genes between IVM-exposed C. elegans N2 and DA1316 strain.  [34]. In addition, eight other DEGs from the DA1316 strain mapped to the QTL (S8 Table). The VL locus contained 32 DEGs (S7 Table) (Fig 6). Six DEGs mapped to the VC locus (Fig 6), which comprised of three stress response genes (F31F7.1, mtl-1, clec-7), a transport-related gene (col-143), a metabolic gene (lipl-4), and an unknown gene (S7 Table). Genes F31F7.1, mtl-1, clec-7, lipl-4, and spp-27 were downregulated between 1.5 and 3.6-fold, with the exception of col-143, which was upregulated 2.4-fold (Fig 6  & S7 Table). Seven DEGs mapped to the VR locus (Fig 6) composed of a GluCl subunit (glc-1),

PLOS ONE
a transporter (F11A5.9) and five unknown genes (S7 Table). The expression of the genes glc-1, F11A5.9, and oac-18 was decreased by 1.4-fold, 1.4-fold, and 2.6-fold, respectively. The remaining genes had 1.54-to 10-fold higher expression levels (Fig 6). A list of candidate genes for further study (S9 Table) has been created, including genes with an opposite expression between N2 and DA1316 C. elegans strain, genes within Abamectin-QTL, exclusive genes in DA1316 strain, and other putative drug-responsive genes.

Discussion
Parasitic nematodes pose a threat to human and animal health [2,3], while the development of resistance against MLs such as IVM exacerbates the problem [16,17,63]. Therefore, it is crucial to understand the molecular mechanisms underlying ML resistance. In this study, we investigated gene expression in wild-type (N2) IVM-exposed C. elegans adults at maximum reproduction to understand the transcriptomic response to IVM exposure. Furthermore, the transcriptomic data was compared to microarray data from the IVM-resistant C. elegans (DA1316) strain [34], and the recently identified Abamectin-QTL [35]. We identified 615 differentially expressed genes in the wild-type C. elegans (N2) strain following exposure to IVM, of which 95% occurred in the IVM 10 −7 M concentration. Only 31 N2 strain DEGs overlapped with the IVM-resistant C. elegans (DA1316) strain, 19 of which, including folate transporter FOLT-2 (folt-2) and transmembrane transporter T22F3.11 (T22F3.11), displayed an opposite expression in either strain. Eighteen of the overlapping genes were enriched for metabolic Phase II enzymes, transmembrane transport, and stress response processes. In addition, we identified 45 differentially expressed genes that mapped within the Abamectin-QTL on C. elegans chromosome V, nine of which overlapped with DEGs from the IVM-resistant C. elegans (DA1316) strain. Following comprehensive analysis we identified potential candidate genes, T-type calcium channel CCA-1 (cca-1), potassium chloride cotransporter KCC2 (kcc-2), folate transporter FOLT-2 (folt-2) and alpha subunit of glutamate-gated channel GluCl-1 (glc-1) for further investigation with possible involvement in IVM resistance.
The ML drug target glc-1 and the putative drug target lgc-26 encoding ion channel subunits were downregulated in the IVM 10 −7 M concentration. The significance of glc-1 in ML response is evident through RNAi knockdown and loss-of-function-mutations [26, 36,64,65]. Heterologous expression of C. elegans glc-1 and formation of a functional chloride channel has been reported in Xenopus oocytes [66]. The role of lgc-26 in ML response remains unknown in C. elegans, but other cys-loop GABA receptor members have been implicated in the ML response of parasitic nematodes. For example, lgc-37 was upregulated in the horse roundworm parasite, Parascaris univalens after IVM exposure [41] and the lgc-54 in the sheep parasite Teladorsagia circumcincta was reported as potential candidate in IVM resistance through genome-wide studies [67]. However, Evans, Wit [35] recently challenged this claim, as they observed that lgc-54 mutants did not display a competitive advantage over wild-type in control conditions.
We observed a downregulation of the calcium channel alpha subunit (cca-1) and potassium chloride cotransporter (kcc-2) genes after IVM exposure. Gene cca-1 regulates pharyngeal pumping by facilitating the effective start of action potentials by influx of calcium ions in response to marginal cell motor neuron stimulation in C. elegans [68,69]. A loss of function in cca-1 has been reported to significantly reduce pharyngeal pumping [69]. Furthermore it is established that IVM reduces pharyngeal pumping in C. elegans by allosterically modulating GluCl hence increased influx of chloride ions [9]. Hypothetically, a pharyngeal pump regulator such as cca-1 would be upregulated in response to decreased pharyngeal pumping rate, however, our results showed the opposite trend. Alternatively, the downregulation of cca-1 could be a downstream effect, a consequence of the negative membrane potential induced by the chloride ion influx, or a defense/ protective mechanism aimed at reducing the adverse effects of IVM. Nonetheless, this highlights the importance of cca-1 in IVM response in C. elegans. The kcc-2 gene encodes a potassium chloride cotransporter, which, in conjunction with the sodium-driven chloride-bicarbonate transporter abts-1, mediates inhibitory GABA signaling. Mediation involves control of cellular chloride gradient to maintain membrane potential in neurons that control locomotion [70]. In addition, kcc-2 and abts-1 mutant of C. elegans were reported to exhibit paralysis after exposure to a GABA receptor agonist, muscimol [70]. Considering that IVM agonistically mediates GABA signaling through the inflow of chloride ions to cause paralysis [9], the role of kcc-2 in IVM response is unknown and warrants exploration. Overall, these findings shed light on additional channels, such as cca-1 and kcc-2, which could be investigated further as potential candidates for involvement in IVM resistance.
Seven of the 11 differentially expressed Transcription factors were downregulated and predominantly enriched for NHRs, which included nhr-17, nhr-21, nhr -55, nhr-57, nhr-99, nhr-115, nhr-119, nhr-121, nhr-137 and nhr-173. Although NHRs have generally been suggested to regulate expression of detoxification genes (phase I, II and III), only a handful have experimental evidence (see review [61]). To our knowledge, none of the NHRs in the current study has been implicated in xenobiotic detoxification regulation. However, a notable example, nhr-8 upregulation has been reported to increase expression of detoxification genes cyp14A2, cyp14A5 and cyp37B1, and Pgp genes (pgp-1,-3,-6,-9 and -13) resulting in the reduced IVM efficacy [71]. Recent work by Guerrero et al. [72] emphasized this phenomenon, demonstrating nhr-8 role in upregulation of Pgp genes (pgp-3, -5, - 11 and -13) in the presence of the drug tunicamycin. The GATA-type transcription factor, elt-2 was downregulated in the current study. The proposed genes regulated by elt-2 in C. elegans include those involved in protection against xenobiotic compounds: CYPs, GSTs, and UGTs such as ugt-22 [73]. While it is conceivable that this downregulation of elt-2 might explain the previously described downregulation of ugt-22, it is unclear why elt-2 is downregulated in the presence of a xenobiotic such as IVM.
In this study, an analysis of putative apoptotic genes showed that four out of the five genes identified were downregulated, including C2H2-type zinc finger transcription factor ces-1 and basic region leucine-zipper (bZIP) transcription factor ces-2. These genes have been previously shown to act as repressors of the key apoptotic activator egl-1 under anti-apoptotic conditions. Under pro-apoptotic conditions, egl-1 represses ced-9, thereby activating the pro-apoptotic genes ced-3 and ced-4, hence apoptosis. The downregulation of ces-1 and ces-2 in this study is indicative of egl-1 being liberated and subsequently able to repress ced-9, thereby activating the pro-apoptotic genes ced-3 and ced-4, leading to apoptosis [62]. Additionally, caudal-type homeodomain transcription factor pal-1, a reported activator of ced-3 [74], was found to be upregulated. To further support the pro-apoptotic response hypothesis, genes ces-1, ces-2, and pal-1 were predominantly present in the highest IVM concentration (10 −7 M), in which all worms were previously reported to be immobile [41]. While these findings suggest that IVM at concentrations of 10 −7 M and 10 −8 M partially activates the apoptotic machinery in the N2 C. elegans strain, it is important to note that other mechanisms involving above genes may also be at play as not all relevant apoptotic genes were detected in the data.
Our second objective was to compare C. elegans strains N2 (wild-type) and DA1316 (IVMresistant) following IVM exposure. Comparing the RNAseq and microarray data from the N2 and DA1316 strains, respectively, revealed that the N2 strain had four times the number of differentially expressed genes with only a 4% (31 genes) overlap between the two strains. Nineteen of these overlapped genes displayed an opposite expression and may have a role in IVM metabolism or resistance. For example, the folate transporter folt-2, whose expression was reduced 6-fold and 2-fold in the resistant DA1316 strain, but increased 4-fold in the sensitive N2 strain, may be involved in the import of IVM. Another potential candidate is the transmembrane transporter T22F3.11, whose expression increased by 3-fold in the resistant strain, but decreased 7-fold in the sensitive N2 strain, suggesting a role in IVM efflux. In addition, genes are that are exclusively expressed in the resistant strain could also be an important group for further exploration. Overall, further studies are needed to confirm the role of these genes in IVM response. Twelve of the 31 overlapping genes were downregulated, including four pud-genes, pud-2.1, pud-3, and pud-4. Pud-genes are unique to Caenorhabditis spp. without known orthologues in other nematodes and their function remains elusive. In spite of that, Cui et al. [75] reported slow growth and hypersensitivity to cadmium of C. elegans after pud-4 (F15E11.12) RNAi knockdown. Another study reported downregulation of pud-1.2 and pud-4 following viral infection of C. elegans [76]. Other overlapping genes such as hsp-17 encoding a heat shock protein, exhibited inverted regulation in either strain. Gene hsp-17, a stress responder, was downregulated 1.7-fold in DA1316 and upregulated 4-fold in N2 strain. This re-emphasizes the inherent sensitivity of the N2 strain to IVM than its resistant or tolerant counterpart. The discrepancy in gene expression between the two studies could be attributable to strain, technique, life-stage, or IVM concentration. Strain differences may account for most of the difference in DEGs, since N2 is more sensitive to IVM than resistant or tolerant strains [34,77]. Thus, the N2 strain may have more elevated cellular and biological processes, which could account for the increase in DEGs. Overall, we speculate that pud-genes may be universal stress responders exclusively in C. elegans.
As a final objective, we mapped the 615 DEGs derived from the IVM-exposed C. elegans to the recently identified Abamectin-QTL [35]. We identified 45 DEGs that corresponded to the QTL regions and compared them to the proposed list of candidate genes from the QTL [35], but only two genes, putative folate transporter (folt-2) and dod-3, between the two data sets overlapped. Gene folt-2 supposedly participates in the active uptake of folate [78], a Bclass vitamin central in the synthesis of nucleotides and amino acids [79]. Like other eukaryotes, C. elegans are intrinsically deficient in folate [80], and therefore rely on dietary sources such as OP50 E. coli. In our data, folt-2 was upregulated, which may be caused by restriction in dietary intake due IVM-induced decreased pharyngeal pumping [9]. This could lead to folate deficiency in the worm. Whether this folate deficiency triggers upregulation of folt-2 remains unclear and therefore entails further study. In contrast, the function dod-3 and its role ML response is unknown. Overall, the role of folt-2, dod-3 and other 43 genes that corresponded to the Abamectin-QTL in ML response is unknown and prompts further investigation.
In conclusion, we have profiled the transcriptomes of adult wild-type (N2) C. elegans after exposure to different IVM concentrations and revealed predominantly downregulated diverse sets of genes. Some genes overlapped with a previous microarray study on IVM-resistant C. elegans while others mapped to recently described abamectin-QTL. Based on this, we have a compiled a list of potential candidate genes for further investigation. Although Abamectin and Ivermectin (IVM) belong to the same subgroup of avermectins, it is important to exercise caution when making comparisons between the two drugs. While Abamectin serves as a precursor drug for IVM, the pharmacokinetics of these drugs have been reported to be considerably different [81,82]. These differences have been attributed, in part, to variations in the drugs' lipophilicity [83]. Future studies employing quantitative genetics approaches, such as identifying IVM-QTL and corresponding candidate genes, validated through the use of near-isogenic lines, mutant strains, and transcriptomic profiling, would provide more detailed information.  Table. Over-representation analysis of IVM-exposed C. elegans (N2) differentially expressed genes.